%% 第6章 信号变换
f=str2sym('sin(a*k)');
F=ztrans(f)
F=str2sym('z/(z-3)^2');
f=iztrans(F)

%% example 6-2
clear,clc,close all;clc;
B=poly([-0.4 -0.4]);
A=poly([-0.8 -0.8 0.5 0.5 -0.1])
[R P K]=residuez(B,A);
R=R',P=P',K

%% 绘制H(z)=(z-1.3)/z的频响曲线
B=[1 -1.3];
A=[1 0];
[H,w] = freqz(B,A,400,'whole');
Hf=abs(H);
Hx=angle(H);
subplot(121);plot(w,Hf);title('离散系统幅频特性曲线');grid on 
xlabel('频率');ylabel('幅度');
subplot(122);plot(w,Hx);title('离散系统相频特性曲线');grid on 
xlabel('频率');ylabel('幅度');
%% example 6-4
clear,clc,clf
w=(-4*pi:0.001:4*pi)+eps;
X=1./(1-0.6*exp(-j*w));
subplot(211);plot(w/pi,abs(X));title('幅频响应');
xlabel('\omega/\pi');ylabel('|H(e^j^\omega)|')
axis([-3.2 3.2 0.5 2.2]);grid
subplot(212);plot(w/pi,angle(X));title('相频响应');
xlabel('\omega/\pi');ylabel('\theta(\omega)')
axis([-3.2 3.2 -0.6 0.6]);grid
set(gcf,'color','w');

%% 离散系统函数零点分析
a=[3 -2 0 0 0 1];
b=[2 1];
ljdt(a,b)
p=roots(a)
q=roots(b)
pa=abs(p)

%% 频率分辨率
% 增加频率分辨率只能通过增加取样点N，但不能通过补零来增加N，而在时域上补零等于在频域上
% 进行插值
clear,clc,clf
Dt=0.00005;t=-0.005:Dt:0.005; % 模拟信号
xa=exp(-2000*abs(t));
Ts=0.0002;n=-25:25;
x=exp(-1000*abs(n*Ts));
K=500;k=0:1:K;w=pi*k/K;
X=x*exp(-j*n'*w);X=real(X);
w=[-fliplr(w),w(2:501)];
X=[fliplr(X),X(2:501)];
% figure
subplot(221);plot(t*1000,xa,'.');title('离散信号');hold on
ylabel('x1(t)');xlabel('t');
stem(n*Ts*1000,x);hold off
legend('x','xa');
subplot(222);plot(w/pi,X,'.');title('离散傅里叶变换');
ylabel('X1(jw)');xlabel('f');
Ts=0.001;n=-25:25;
% 离散时间信号
x=exp(-1000*abs(n*Ts));
K=500;k=0:1:K;w=pi*k/K;
% 离散时间傅里叶变换
X=x*exp(-j*n'*w);X=real(X);
w=[-fliplr(w),w(2:501)];
X=[fliplr(X),X(2:501)];
subplot(223);plot(t*1000,xa,'.');title('离散信号')
xlabel('t');ylabel('x2(t)');hold on
stem(n*Ts*1000,x);hold off
subplot(224);plot(w/pi,X,'.');title('离散时间傅里叶变换')
xlabel('f');ylabel('X2(jw)');

%% 6-13
clear,clc
fs=100;     %采样频率
Ndata=32;   %数据长度
N=32;       %FFT数据长度
n=0:Ndata-1;
t=n/fs;     %数据对应的时间序列
x=0.5*sin(2*pi*20*t)+2*sin(2*pi*60*t);  %时间域信号
plot(t,x)
y=fft(x,N);  %信号傅里叶变换
mag=abs(y); %求取振幅
f=(0:N-1)*fs/N; %真实频率
subplot(221);plot(f(1:N/2),mag(1:N/2)*2/N);
xlabel('频率/Hz');ylabel('振幅');
title('Ndata=32;FFT所以采点个数=32');grid on

Ndata=32;   %数据长度
N=128;       %FFT数据长度
n=0:Ndata-1;
t=n/fs;     %数据对应的时间序列
x=0.5*sin(2*pi*20*t)+2*sin(2*pi*60*t);  %时间域信号
y=fft(x,N);  %信号傅里叶变换
mag=abs(y); %求取振幅
f=(0:N-1)*fs/N; %真实频率
subplot(222);plot(f(1:N/2),mag(1:N/2)*2/N);
xlabel('频率/Hz');ylabel('振幅');
title('Ndata=32;FFT所以采点个数=128');grid on

Ndata=136;   %数据长度
N=128;       %FFT数据长度
n=0:Ndata-1;
t=n/fs;     %数据对应的时间序列
x=0.5*sin(2*pi*20*t)+2*sin(2*pi*60*t);  %时间域信号
y=fft(x,N);  %信号傅里叶变换
mag=abs(y); %求取振幅
f=(0:N-1)*fs/N; %真实频率
subplot(223);plot(f(1:N/2),mag(1:N/2)*2/N);
xlabel('频率/Hz');ylabel('振幅');
title('Ndata=136;FFT所以采点个数=128');grid on

Ndata=136;   %数据长度
N=512;       %FFT数据长度
n=0:Ndata-1;
t=n/fs;     %数据对应的时间序列
x=0.5*sin(2*pi*20*t)+2*sin(2*pi*60*t);  %时间域信号
y=fft(x,N);  %信号傅里叶变换
mag=abs(y); %求取振幅
f=(0:N-1)*fs/N; %真实频率
subplot(224);plot(f(1:N/2),mag(1:N/2)*2/N);
xlabel('频率/Hz');ylabel('振幅');
title('Ndata=136;FFT所以采点个数=512');grid on

%% x(n)是一个正弦信号、一个余弦信号和白噪声叠加的信号，试用FFT对其进行频谱分析
clc,clear
N=1024;f1=0.1;f2=0.2;fs=1;
a1=5;a2=3;
w=2*pi/fs;
x=a1*sin(w*f1*(0:N-1))+a2*cos(w*f2*(0:N-1))+randn(1,N); %信号
% 应用FFT求频谱
subplot(211);plot(x(1:N/4));title('原始信号');
f=-0.5:1/N:0.5-1/N;
X=fft(x);
subplot(212);plot(f,abs(X));title('频域信号');
%% 逆FFT
clear,clc
fs=200;
N=128;
n=0:N-1;
t=n/fs;
x=0.5*sin(2*pi*20*t)+2*sin(2*pi*60*t);
subplot(221);plot(t,x);title('原始信号');grid on
xlabel('时间/s');ylabel('x');
y=fft(x,N);
mag=abs(y);
f=n*fs/N;   %频率序列
subplot(222);plot(f(1:N/2),mag(1:N/2)*2/N);title('FFT');
xlabel('频率/Hz');ylabel('振幅');
xifft=ifft(y);
realx=real(xifft);
ti=[0:length(xifft)-1]/fs;      %时间序列
subplot(223);plot(ti,realx);title('IFFT后的信号');grid on
xlabel('时间/s');ylabel('x');
yif=fft(xifft,N);
mag=abs(yif);
f=[0:length(y)-1]'*fs/length(y);
subplot(224);plot(f(1:N/2),mag(1:N/2)*2/N);title('IFFT后的快速傅里叶变换');
xlabel('频率/Hz');ylabel('振幅');

%% 利用FFT对叠加白噪声信号进行滤波，将4~8Hz的波滤除
clear,clc,clf
dt=0.05;N=1024;
n=0:N-1;t=n*dt;
f=n/(N*dt); %频率序列
f1=3;f2=10;%两个正弦信号的频率
x=0.5*cos(2*pi*f1*t)+sin(2*pi*f2*t)+randn(1,N);
subplot(2,2,1);plot(t,x)
title('原始信号时间域');xlabel('t/s');
y=fft(x);
xlim([0 12]);ylim([-2 2]);
subplot(2,2,2);plot(f,abs(y)*2/N);
title('原始振幅谱');
xlabel('频率/Hz');ylabel('振幅');
xlim([0 25]);ylim([0 0.6]);

f1=4;f2=8;      %要滤除的频率上下限
yy=zeros(1,length(y));  %设置与y相同的元素数组
for m=0:N-1 
    if(m/(N*dt)>f1 && m/(N*dt)<f2) %|| (m/(N*dt)>(1/dt-f2) && m/(N*dt)<(1/dt-f1))
        yy(m+1)=0;
    else
        yy(m+1)=y(m+1);
    end
end
subplot(2,2,3);plot(t,real(ifft(yy)));
xlabel('t/s');ylim([-1.5 1.5]);xlim([0 12]);
subplot(2,2,4);plot(f,abs(yy)*2/N)
xlim([0 25]);ylim([0 0.6]);
xlabel('频率/Hz');ylabel('振幅');
title('滤除后的振幅谱');
%% 对信号进行离散余弦变换
clear,clc,clf
n=1:100;
x=10*sin(2*pi*n/20)+20*cos(2*pi*n/30);
y=dct(x);
subplot(1,2,1);plot(x);title('原始信号');
subplot(1,2,2);plot(y);title('DCT信号');




